#!/usr/bin/perl -w
#-----------------------------------------------------------+
#                                                           |
# TARGET SITE DUPLICATION SELECTION                         |
#                                                           |
#-----------------------------------------------------------+
#  AUTHOR: James C. Estill                                  |
# CONTACT: jestill_at_sourceforge.net                       |
# STARTED: 06/08/2006                                       |
# UPDATED: 06/08/2006                                       |
#                                                           |
# DESCRIPTION:                                              |
#  Extract target site duplication sequences from a genomic |
#  sequence given the start and stop positions of an LTR    |
#  retrotransposon.                                         |
#                                                           |
# INPUT VARIABLES:                                          |
#   [f] FASTA file containing ALL contigs of interest       |
#       This can be a single chrom, multiple chrom,         |
#       multiple BACs etc.                                  | 
#   [a] Annotation file in GFF Format                       |
#   [o] Output file path                                    |
#   [l] Length of LTR seq to return                         |
#       Default value is 3                                  |
#   [t] Length of putative TDS to return                    |
#       Default value is 5                                  |
#   [q] Quiet mode. Do not print program status to screen   | 
#                                                           |
# REQUIRES:                                                 |
#   - bioperl                                               |
#     Makes use of the BioPERL seqIO function               |
#                                                           |
# USAGE:                                                    |
#  tsd_sel.pl -f SeqFile.fasta -a LTRAnn.gff -o OutFile.txt |
#-----------------------------------------------------------+

#-----------------------------+
# INCLUDES                    |
#-----------------------------+
use Bio::SeqIO;                # Allows for fasta sequence input/output
use Term::ANSIColor;           # Allows for pretty print of ANSI output
use Getopt::Std;               # Allow for option input from cmd line

#-----------------------------+
# USER SET VARIABLES          |
#-----------------------------+
# The following set will be used temporarily to get the 
# basic functions of the program working. These then
# be changed to cmd line variables. The choice of random
# genes or a select set of genes will left up to the
# user.
my $WorkDir = "/home/jestill/projects/clem/";
my $LogFile = $WorkDir."TSDNotes.txt";
my $sep="|";                   # Character to separate the feature name from
                               # the clone name in the fasta output file.
                               # Feature name is name of the LTR.

#-----------------------------+
# GLOBAL VARIABLES            |
#-----------------------------+
# These variables need to be local to the program
# but need to be globally avaiable to the 
# subfunctions. The user will not need to 
# modify these variables.
my @StartGeneExons;            # Array to hold information for start exon
my @Annot;                     # Array to hold all annotation info for BAC
my @StartGeneList;             # Array to hold the list of gene names
                               # to use as the start genes.
my $StartPos;                  # Random start position on the BAC for
                               # one end of a simulated clone
my $SearchPos;                 # Search position to look for gene on other
                               # end of the fake clone.
my $CumSum = "0";              # Var to hold cumulative sum
my $ans;                       # Var to hold answers to questions
my $NumExonStartGene;          # Number of exons in the gene model
my $UsableLen;                 # Usable length of the exon
my $SeqId;                     # Sequence ID of contig sequence record to use
my $RandCumSum;                # Random cumulative sum
my $SeqLength;                 # Length of the sequence file
my $NumTSDMatch;               # Number of matches in TSD
my $MatchString;               # String to visualize match
my $MaxNum;                    # Max position of LTR elements when
                               # stored in the array

#-----------------------------+
# GET VARIABLES FROM CMD LINE |
#-----------------------------+
# Declare the Options hash to hold the options from the command line
my %Options;
getopts('f:a:o:t:l:q', \%Options);

my $FastaFile = $Options{f} || 
    die "\nERROR: A fasta file was not listed with the -f flag\n";
my $AnnotFile = $Options{a} ||
    die "\nERROR: An annotation file was not listed with the -a flag\n";
my $OutFile = $Options{o} ||
    die "\nERROR: An output file was not listed with the -o flag.\n";

my $TSDLen = $Options{t} || "5";
my $LTRLen = $Options{l} || "3";
my $quiet = $Options{q};


#-----------------------------+
# LET THE USER SEE THE VARS   |
# AND CHOOSE TO CONTINUE TO   |
# EXIT                        |
#-----------------------------+
&UserVarCheck;

#-----------------------------+
# FILEHANDLES I/O             |
#-----------------------------+
open (OUTFILE, ">$OutFile") ||
    die ("Can not open outfile:\n$OutFile");

open (LOGFILE, ">$LogFile") ||
    die ("Can not open logfile:\n$LogFile");

my $seqio_object = Bio::SeqIO->new( '-file' => $FastaFile,
				    '-format' => "fasta" ) ||
    die ("Can not open fasta file:\n$FastaFile");


#-----------------------------+
# LOAD ANNOTATION DATA        |
#-----------------------------+
&LoadAnnot($AnnotFile);        # Load full annotation file to 2d Array
my $LenAnnot = @Annot; 
print "$LenAnnot Annotation Objects\n";

#-----------------------------+
# FOR EACH SEQUENCE RECORD    |
# IN THE FASTA FILE SEARCH FOR|
# LTRS TO EXTRACT             |
#-----------------------------+
while (my $inseq = $seqio_object->next_seq)
{

    $SeqId = $inseq->primary_id;
    print color 'bold';
    print "SEQUENCE: ".$SeqId."\n";
    print color 'reset';

    #-----------------------------+
    # FOR EACH ANNOTATED LTR      |
    # IN THE GFF FILE             |
    #-----------------------------+
    $MaxNum = $LenAnnot - 1;
    for ($i=0; $i<=$MaxNum; $i++)
    {
	if ($SeqId =~ $Annot[$i][0])
	{
	    my $LTR_Seq = $Annot[$i][0];
	    my $LTR_ID = $Annot[$i][8];
	    my $LTR_Start = int( $Annot[$i][3] );
	    my $LTR_End = int( $Annot[$i][4] );
	    print  "\tLTR:".$LTR_ID." : ".
		$LTR_Start."-".$LTR_End."\n";

	    my $LeftEnd = $inseq->subseq( $LTR_Start - $TSDLen, 
					  $LTR_Start + $LTRLen-1);

	    my $RightEnd = $inseq->subseq( $LTR_End - $LTRLen+1,
					   $LTR_End + $TSDLen );

	    my $LeftMatch = substr($LeftEnd, 0, $TSDLen);
	    my $RightMatch = substr($RightEnd, $LTRLen);
	    my $MatchLen = length($LeftMatch);

	    #-----------------------------+
	    # GENERATE A STRING SHOWING   |
	    # MATCHING RESIDUES           |
	    #-----------------------------+
	    $NumTSDMatch = "0";
	    $MatchString = "";
	    for ($j=0; $j<=$MatchLen; $j++)
	    {

		my $L = substr($LeftMatch,$j,1);
		my $R = substr($RightMatch,$j,1);
		
		if ($L =~ $R)
		{
		    $NumTSDMatch++;
		    $MatchString = $MatchString."|";  
		}else{
		    $MatchString = $MatchString." "; 
		}

	    }
	    
	    #-----------------------------+
	    # PRETTY PRINT OUTPUT TO THE  |
	    # SCREEN                      |
	    # THIS WILL BE SORTED BY THE  |
	    # BAC OCCURRENCE IN THE FASTA |
	    # INPUT FILE                  |
	    #-----------------------------+
	    if ($LeftMatch =~ $RightMatch)
	    {
		print color 'green';
		print "\t\tMATCH:".$NumTSDMatch."\n";
		print "\t\t";
		print " " x $LTRLen;
		print $LeftEnd."\n";
		print "\t\t";
		print " " x $LTRLen;
		print $MatchString."\n";
		print "\t\t".$RightEnd."\n";
		print color 'reset';
	    }else{
		print color 'red';
		print "\t\tMATCH:".$NumTSDMatch."\n";
		print "\t\t";
		print " " x $LTRLen;
		print $LeftEnd."\n";
		print "\t\t";
		print " " x $LTRLen;
		print $MatchString."\n";
		print "\t\t".$RightEnd."\n";
		print color 'reset';
	    }


	    #-----------------------------+
	    # SAVE RESULTS TO ANNOTATION  |
	    # ARRAY                       |
	    #-----------------------------+
	    $Annot[$i][9] = $NumTSDMatch;
	    $Annot[$i][10] = $LeftEnd;
	    $Annot[$i][11] = $RightEnd;
	    $Annot[$i][12] = $MatchString;
	    $Annot[$i][13] = $LeftMatch;
	    $Annot[$i][14] = $RightMatch;

	} # End of found seq/annotation match

    } # End of for each annotation in the gff file

} # End of for each sequence in the output file


#-----------------------------+
# PRINT ANNOTATION RESULTS    |
# TO TAB DELIM OUTPUT FILE    |
# THIS WILL BE IN THE         |
# ORIGINAL INPUT ORDER        |
#-----------------------------+
print "\n\n";
print "==================================\n";
print "   NOW SHOWING RESULTS IN ARRAY\n";
print "      IN A PLAIN PRINT FORMAT\n";
print "==================================\n";

print color 'bold';
print "BAC\tELEMENT\tMATCH\tLEFT\tRIGHT\t\n";
print color 'reset';

for ($i=0; $i<=$MaxNum; $i++)
{
    print $Annot[$i][0]."\t";
    print $Annot[$i][8]."\t";
    print $Annot[$i][9]."\t";
    #print $Annot[$i][10]."\t";
    #print $Annot[$i][11]."\t";
    print $Annot[$i][13]."\t";
    print $Annot[$i][14]."\t";
    print "\n";
    
}

#-----------------------------+
# PRINT RESULTS IN THE        |
# ORIGINAL ORDER TO THE       |
# OUTFILE                     |
#-----------------------------+
for ($i=0; $i<=$MaxNum; $i++)
{
    print OUTFILE $Annot[$i][8]."\t";
    print OUTFILE $Annot[$i][0]."\t";
    print OUTFILE $Annot[$i][9]."\n";
    print OUTFILE "\t\t";
    print OUTFILE " " x $LTRLen;
    print OUTFILE $Annot[$i][10]."\n";
    print OUTFILE "\t\t";
    print OUTFILE " " x $LTRLen;
    print OUTFILE $Annot[$i][12]."\n";
    print OUTFILE "\t\t".$Annot[$i][11]."\n";
    print OUTFILE "\n";    
}


close (OUTFILE);
close (LOGFILE);

exit;

#-----------------------------------------------------------+
# SUBFUNCTIONS                                              |
#-----------------------------------------------------------+

#-----------------------------+
# LOAD BAC ANNOTATION DATA    |
#-----------------------------+
# This should be in tab delim gff format.
sub LoadAnnot 
{
    my $GffFile = $_[0];
    my (@line, $SeqName, $Source, $Feature, $Start, $End);
    my ($Score, $Strand, $Frame);

    open (INFILE, $GffFile) || 
	die ("Can not load file $GffFile \n");
    while (<INFILE>)                # For each line in the file
    {

	# Get in from tab delim text file 
	my @line = split(/\t/, $_); 
	my $Seqname = $line[0];
	my $Source = $line[1];
	my $Feature = $line[2];
	my $Start = $line[3];
	my $End = $line[4];
	my $Score = $line[5];
	my $Strand = $line[6];
	my $Frame = $line[7];
	my $Attribute = $line[8]; # Gene name column
	chomp $Attribute;         # Get rid of newline character

	# Can get rid of comment below to show the loading
	# of the individual features
	#print "Load: $Attribute $Start\-$End\n";

	# Load the information into the @Annot array
	push (@Annot, [$Seqname, $Source, $Feature, $Start, $End,
		       $Score, $Strand, $Frame, $Attribute]);

    }
    close (INFILE);

    my $TestAryLen = @Annot;
    if ($TestAryLen < 1)
    {
	print color 'red';
	print "ERROR: I did not find any annotation data\n";
	print color 'reset';
	exit;
    }

}

#-----------------------------+
# USER FEEDBACK SUBFUNCITON   |
#-----------------------------+
sub UserFeedback
{
  
  my $Question = $_[0];
  my $Answer;
  
  print "\n$Question \n";
  
  while (<>)
    {
      chop;      # Get rid of the carriage return character
      if ( ($_ eq 'y') || ($_ eq 'Y') || ($_ eq 'yes') || ($_ eq 'YES') )
	{
	  $Answer = "Y";
	  return $Answer;
	}
      elsif ( ($_ eq 'n') || ($_ eq 'N') || ($_ eq 'NO') || ($_ eq 'no') ) 
	{
	  $Answer = "N";
	  return $Answer;
	}
      else
	{
	  print "\n$Question \n";
	}
    }

}

#-----------------------------+
# USER VARIABLE CHECK         |
#-----------------------------+
# Let's the user check the input variables

sub UserVarCheck
{
    print "\nYOU HAVE SELECTED THE FOLLOWING VARIBLES:\n";
    #-----------------------------+
    # FASTA FILE                  |
    #-----------------------------+
    print "FASTA File:\n";
    if (-e $FastaFile)
    {
	print color 'green';
	print "\t$FastaFile\n";
	print color 'reset';
    }else{
	print color 'red';
	print "\tWARNING: $FastaFile \n\tDoes Not Exist\n";
	print color 'reset';
    }
    
    #-----------------------------+
    # ANNOTATION FILE             |
    #-----------------------------+
    print "ANNOTATION FILE:\n";
    if (-e $AnnotFile)
    {
	print color 'green';
	print "\t$AnnotFile\n";
	print color 'reset';
    }else{
	print color 'red';
	print "\tWARNING: $AnnotFile \n\tDoes Not Exist\n";
	print color 'reset';
    }
    
    #-----------------------------+
    # OUTPUT FILE                 |
    #-----------------------------+
    print "OUTPUT FILE:\n";
    if (-e $OutFile)
    {
	print color 'red';
	print "\t$OutFile already exists\n";
	print "\tThe existing file will be overwritten\n";
	print color 'reset';
    }else{
	print color 'green';
	print "\t$OutFile\n";
	print color 'reset';
    }
    
    #-----------------------------+
    # LTR LENGTH                  |
    #-----------------------------+
    print "LTR LENGTH:\n\t$LTRLen\n";

    #-----------------------------+
    # TSD LENGTH                  |
    #-----------------------------+
    print "TSD LENGTH\n\t$TSDLen\n";
    
    #-----------------------------+
    # QUIET MODE                  |
    #-----------------------------+
    # If this is in quiet mode the answer to
    # the question is Y and the user will not
    # need to supply feedback at the command line.
    if (! $quiet){
	$ans = &UserFeedback("\nDo you want to continue (y/n)?");
    }else{
	$ans = "Y";
    }
    
    if ($ans =~"N" )
    { 
	print "Goodbye\n";
	exit;
    }elsif ($ans =~"Y"){
	print "Starting the process...\n";	
    }    
}

#-----------------------------------------------------------+
# HISTORY                                                   |
#-----------------------------------------------------------+
#
# 06/08/2006
# - Begin modification of DblHit to do a dts_sel for clem
#
#-----------------------------------------------------------+
# MODEL                                                     |
#-----------------------------------------------------------+
#
# Given a genomic sequence and Start/End Positions of the
# LTR Retrotransposon for a given genomic sequence:
#
# GENOME
# |=========================================================================|
#                    
#                  START                           END
#                   ||                             ||
# LTR               \/                             \/
#                    [=5'LTR=]-------------[3'=LTR=]
#                                      
# RETURNED          
# SEQS              =:                              :=
#                  5'Seq                             3'Seq
#
# GIVEN: 
# START is position of start of 5' LTR
# END is position of end of 3' LTR
